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Abstract: The estimation of small probabilities of failure from computer simulations is a classical problem 
in engineering, and the Subset Simulation algorithm proposed by Au & Beck (Prob. Eng. Mech., 2001) has 
CN ■ become one of the most popular method to solve it. Subset simulation has been shown to provide significant 
savings in the number of simulations to achieve a given accuracy of estimation, with respect to many other 
CN ■ Monte Carlo approaches. The number of simulations remains still quite high however, and this method can be 
^ . impractical for applications where an expensive-to-evaluate computer model is involved. 

. We propose a new algorithm, called Bayesian Subset Simulation, that takes the best from the Subset Simulation 
algorithm and from sequential Bayesian methods based on kriging (also known as Gaussian process modeling). 
Q | The performance of this new algorithm is illustrated using a test case from the literature. We are able to report 
; promising results. In addition, we provide a numerical study of the statistical properties of the estimator. 
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> 1 1. INTRODUCTION 

m 

■ In this paper, we propose an algorithm called Bayesian Subset Simulation (BSS), that combines the Bayesian 
. decision-theoretic framework from our previous studies [1,2] with the Subset Simulation algorithm [3]. 

o: 

£N| . Let r = {x G X : f{x) > u} denote the excursion set of a function / : X — > R above a threshold 
u € R. We are interested in estimating the probability a{u) := Px(r), which corresponds to the probability 
of failure of a system for which / is a function of performance (see, e.g., [2]). If the probability a{u) is 
small, estimating it using the Monte Carlo estimator <5^ c = 1/m YliLi ^ L f(x i )>u> '~ d Px> requires a large 
number of evaluations of /. If the performance function / is expensive to evaluate, this leads to use a large 
amount of computational resources, and in some cases, it may be even impossible to proceed in reasonable 
time. Estimating small probabilities of failure with moderate computational resources is a challenging topic. 

When a(u) is small, the main problem with the estimator a^ c is that the sample size m must be large in 
order to get a reasonably high probability of observing at least a few samples in T. In the literature, importance 
sampling methods have been considered to generate more samples in the failure region Y. However, the success 
of this kind of methods relies greatly on prior knowledge about the failure region T and on a relevant choice for 
the proposal sampling distribution. 

The idea of Subset Simulation is to decompose the difficult problem of generating samples in the failure region 
into a series of easier problems, by introducing intermediate failure events. Let no = — oo < u\ < U2 < ■ ■ ■ < 
ut = u be a sequence of increasing thresholds and define a corresponding sequence of decreasing excursion 
sets T := X D Ti D • • • D T T := T, where T t := {x G X : f(x) > u t }, t = 1,...,T. Notice that 
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Ft = Hi=i ^j- Then, using the properties 



a(ut+i) = a(itt) Px(r t+ i|r t ) , i>0, 



a(it) can be rewritten as a product of conditional probabilities: 



r-i 



a(«) = P x (r T ) = J] p x(r m |r t ). (2) 

t=o 

Thus, the idea of Subset Simulation is to replace the problem of estimating the small probability a(u) by that 
of estimating the higher conditional probabilities Px(IV|-i|I\), < t < T. 

In [3], a standard Monte Carlo simulation method is used to estimate Px(Ti) = Px(Fi|Fo). For the other 
conditional probabilities, a Markov Chain Monte Carlo method is used to simulate samples in T t , and then 
Px(rVfi|I\) is estimated using a Monte Carlo method. Due to the direct use of Monte Carlo method, the 
number of evaluations needed remains still quite high. For many practical applications where the performance 
function corresponds to an expensive-to-evaluate computer model, it is not applicable. Note that the Subset 
Simulation algorithm has recently caught the attention of the Sequential Monte Carlo (SMC) community: using 
standard tools from the SMC literature, [4] derives several theoretical results about some very close versions of 
the Subset Sampling algorithm. 

In this work, we propose an algorithm that takes advantage of a Gaussian process prior about / in order to 
decrease the number of evaluations needed to estimate the conditional probabilities Px(IV|_i|I\). The Gaussian 
process model makes it possible to assess the uncertainty about the position of the intermediate excursion 
sets Tt, given a set of past evaluation results. The idea has its roots in the field of design and analysis of 
computer experiments (see, e.g., [5, 6, 7, 8, 9, 10, 11]). More specifically, kriging-based sequential strategies 
for the estimation of a probability of failure are closely related to the field of Bayesian global optimization 
[12, 13, 14, 15, 16, 17]. 

The paper is organized as follows. In Section 2, we give a detailed presentation of our new Bayesian Subset 
Simulation algorithm. In Section 3, we apply the algorithm on an example from the literature, and we perform 
numerical simulations to investigate the performance of the proposed algorithm. A comparison with Subset 
Simulation and classical Monte Carlo methods is provided. Finally, we conclude in Section 4. 

Remark. An algorithm involving kriging-based adaptive sampling and Subset Simulation has been recently 
proposed by V. Dubourg and co-authors [18, 19] to address the problem of Reliability-Based Design Optimiza- 
tion (RBDO). Their approach is different from the one proposed in this paper, which addresses the problem of 
reliability analysis. 



2. BAYESIAN SUBSET SIMULATION ALGORITHM 



2.1. Algorithm 

Our objective is to build an estimator of a(ur) from the evaluations results of / at a number of points 
X\, X2, ■ ■ ■ , Xn £ X. Let £ be a random process modeling our prior knowledge about /, and for each n > 0, 
denote by T n the cr-algebra generated by X\, . . . , X n , £(X n ). A natural Bayesian estimator of a{ut) 

using n t evaluations is the posterior mean 

a t = E nt (a(u t )) = E nt l^ >Ut dP x ^j = ^c/tdPx, (3) 
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where g t : x 6 X i— > P nt ( £(x) > ut ) and E n (resp. P„) denotes the conditional expectation (resp. conditional 
probability) with respect to T n (see [2]). Note that, gt(x) can be readily computed for any x using kriging (see, 
e.g., [2]). 

Assume now that Px has a probability density function px and consider the sequence of probability density 
functions q t , < t < T, defined by 

1 

qt(x) = —px(x)gt{x). (4) 
a t 

We can write a recurrence relation similar to (1) for the sequence of Bayesian estimators a t : 

a t = [ g t (x)px(x)dx = a t -i [ 9t ^ q t _ 1 (x)dx. (5) 
J J 9t-i{x) 

The idea of our new algorithm, that we call Bayesian Subset Simulation, is to construct recursively a Monte 
Carlo approximation ax of the Bayesian estimator at, using (5) and sequential Monte Carlo simulation (SMC) 
(see, e.g., [20]) for the evaluation of the integral with respect to q t -\ on the right-hand side. More precisely, 
denoting by m the size of the Monte Carlo sample, we will use the recurrence relation 



<»t = at-i x - ^ - l<t<T, (6) 



where variables . . . , Y™ x are distributed according to 1 the density qt-i, which leads to 

t=0 i=l yty t ' 

The connection between the proposed algorithm and the original Subset Simulation algorithm is clear from the 
similarity between the recurrence relations (1) and (5), and the use of SMC simulation in both algorithms to 
construct recursively a "product-type" estimator of the probability of failure (see also in [20], Section 3.2.1, 
where this type of estimator is mentioned in a very general SMC framework). 

Our choice for the sequence of densities qi,...,qx also relates to the original Subset Simulation algorithm. 
Indeed, note that qt(x) oc E nt (l^> Ut Pjl), and recall that q t oc l^> Ui Px is the distribution used in the Subset 
Simulation algorithm at stage t. (This choice of instrumental density is also used by [22, 23] in the context of 
a two-stage kriging-based adaptive importance sampling algorithm. This is indeed a quite natural choice, since 
qr oc l^> u px is the optimal instrumental density for the estimation of a{u) by importance sampling.) 



2.2. Implementation 



This section gives implementation details for our Bayesian Subset Simulation algorithm, the principle of which 
has been described in the Section 2.1. The pseudo-code for the algorithm is presented in Table 1. 

The initial Monte Carlo sample Yo = {Yq, . . . , Y™} is a set of independent random variables drawn from the 
density qo = px — in other words, we start with a classical Monte Carlo simulation step. At each subsequent 
stage t > 1, a new sample Y t is produced from the previous one using the basic reweight/resample/move 
steps of SMC simulation (see [20] and the references therein). In this article, resampling is earned out using a 
multinomial sampling scheme, and the move step relies on a fixed-scan Metropolis-within-Gibbs algorithm as 
in [3] with a Gaussian-random-walk proposal distribution for each coordinate (for more information on these 
techniques, see, e.g., [24]). 

'By "distributed according to", it is not meant that Y^_ 1 , . . . , Y™i are independent and identically distributed. This is never the 
case in sequential Monte-Carlo techniques. What we mean is that the sample Y t _ 1 , . . . , Y™\ is targetting the density qt-i in the sense 
of, e.g., [21]. 
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Table 1: Algorithm of Bayesian Subset Simulation 



a) Initialize (Stage 0): 

1. Generate a MC sample Y t = {Yq, . . . , Y™}, drawn according to the distribution Px 

2. Initial DoEl n = {(X 1 ,f(X 1 )),..., (X na , f(X na ))} (maximin) 

3. Choose kriging model, estimate parameters kg 

b) At each stage t(t = 1 . . . T): 

1. Compute the kriging predictor /„ and choose threshold u l ~ l 

2. Select and evaluate N t new points using a SUR sampling criterion for the threshold u t . 

3. Update Z n , adjust intermediate threshold Ut-i according to 

4. Regenerate a new sample Y t : 

4.1 reweight: calculate weights: w\ oc gt (Xi-i) / 9t-\{Xt-\) 

4.2 resample: generate a sample j according to weights 

4.3 move: for each i < m, Y t l ^ K (Yj-i^, • ) 

c) The final probability of failure is calculated by 



A number of evaluations of the performance function is done at each stage of the algorithm. This number 
is meant to be much smaller than the size m of the Monte Carlo sample — which would be the number of 
evaluations in the classical Subset Sampling algorithm. For the initialization stage (t = 0), we choose a 
space filling set of points Yo as usual in the design of computer experiments [25]. At each subsequent stage, 
we use Nt iteration of a SUR sampling strategy [2] targeting the threshold u t to select the evaluation points. 
Adaptive techniques to choose the sequence of thresholds and the number of points per stage are presented in 
the following sections. 

Remark 1. The resampling step could most certainly benefit from more elaborate schemes, such as the residual 
resampling scheme [26, 27, 28]. The comparison of resampling schemes is left for future work. 



2.3. Adaptive choice of the thresholds u t 

It can be proved that, for an idealized 2 Subset Simulation algorithm with fixed thresholds uq < u\ < ■ ■ ■ < 
ut = u, it is optimal to choose the thresholds to make all conditional probabilities Px(rt + i|r f ) equal to some 
constant value (see [4], Section 2.4). This leads to the idea of choosing the thresholds adaptively in such a way 
that, in the product estimate 

T m 

4 ubSamp = II^EMitO. 

t=l 1=1 

assuming that Yt , Y t m are independent and identically distributed according to qt. 
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each term but the last is equal to some predefined constant p$. In other words, u t is chosen as the (1 — 
Po)-quantile of Yt-\. This idea was first suggested by [3] in Section 5.2, on the heuristic ground that the 
algorithm should perform well when the conditional probabilities are neither too small (otherwise they are hard 
to estimate) nor too large (otherwise a large number of stages is required). The asymptotic behavior of the 
resulting algorithm, when m is large, has been analyzed by [4]. 

In Bayesian Subset Simulation, we propose to choose the thresholds adaptively using a similar approach. More 
precisely, considering the product form of the estimator (7), we suggest to choose u t in such a way that 

9t{Yi) p °- 

The equation can be easily solved since the left-hand side is a strictly decreasing function of ut- 

Remark 2. Note that [4] proved that choosing adaptive levels in Subset Simulation introduces a positive bias 
of order 1/m, which is negligible compared to its standard deviation. 

2.4. Adaptive choice of the number N t of evaluation at each stage 

In this section, we propose a technique to choose adaptively the number Nt of evaluations of the performance 
function that must be done at each stage of the algorithm. 

Let us assume that t > 1 is the current stage number; at the beginning of the stage, n t -\ evaluations of the 
performance function are available from previous stages. After several additional evaluations, the number of 
available observations of / is n > n t -\. Then, for each i £ {1, . . . , m}, the probability of misclassification 3 
of x E X with respect to the threshold ut is 

T tj n(x) = mm(p n (x,u t ), l-p n (x,,u t j), 

where p n (x, u) = E n (1«j >b ) (see [2]). We shall decide to stop adding new evaluations at stage t when 

i=l 

for some prescribed r\ > 0. 

3. NUMERICAL RESULTS 

In this section, we apply the proposed algorithm on a simple 2D test case from the structural reliability literature. 
The problem under consideration is the deviation of a cantilever beam, with a rectangular cross-section, and 
subjected to a uniform load [29, 30]. The performance function is: 

f(x 1 ,x 2 ) = 18.46154 - 7.476923 x 10 10 ^ . (8) 

X2 A 

The uncertain factors are x\ and X2, which are supposed to be independent and normally distributed, as specified 
in Table 2. We use u = 17.8 as the threshold for the definition of the failure event. The probability of failure, 
which will be used as reference estimator, obtained using S^ c with m = 10 8 , is approximately 3.85 x 10 -5 
(with a coefficient variance of about 1/ yjma « 1.61%). Figure 1 shows the distribution of the input factors 
along with a contour plot of /. Notice that the failure region is quite far from the center region of the input 
distribution. 
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Table 2: Random input factors 



Variable 


Distribution 


Mean m 


Standard deviation a 


Xi 


N 


0.001 


0.0002 


X2 


M 


250 


37.5 




In the Bayesian Subset Simulation algorithm, we set an initial design of size No = 10 which is equal to five 
times the dimension d of the input space (In the literature, very little is known about the problem of choosing 
Nq, however some authors recommend to start with a sample size proportional to the dimension d, see [31]). 
Concerning the choice of Nq, we decide to apply a greedy MAXMIN algorithm and sequentially choose the 
points which will maximize the minimal Euclidean distance between any two points in the initial Monte Carlo 
sample Y(j. At each stage, we choose a Monte Carlo sample of size m = 1000. A Gaussian process with 
constant unknown mean and a Matern covariance function is used as our prior information about /. The 
parameters of the Matern covariance functions are estimated on the initial design by REML (see, e.g., [32]). 
In this experiment, we follow the common practice of re-estimating the parameters of the covariance function 
during the sequential strategy, and update the covariance function after each evaluation. The target conditional 
probability between successive thesholds is set to po = 0.1. The intermediate threshold u t is chosen by the 
criterion (2.3). The proposal distribution q t for a Gibbs sampling is a Gaussian distribution M(0,a 2 ), where 
a 2 is specified in Table 2. The stopping criterion for the adaptive SUR strategy is set to rj t = 10~ 6 (for 
t = 1, . . . , T - 1) and r) T = 10~ 7 . 

Figure 2 shows the Design of Experiment (DoE) selected by the algorithm at stage t = 1, 2, 3 and the last stage 
for one run. Table 3 lists the number of evaluations (rounded to integer) at each stage averaged over 50 runs. 
We can see that an average total of evaluations N = YlJ=o — 104 are needed for our proposed Bayesian 
Subset Simulation, while for Subset Simulation, the number is 1000 + 900 x 4 = 4600. 

3 See [2] Section 2.4 for more information 
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Table 3: Average number of evaluations at each stage. 



N t 


1 


2 


3 


4 


5 


Sub-Sim 


1000 


900 


900 


900 


900 


Bayesian Sub-Sim 


14 


17 


17 


18 


28 




Figure 2: Evaluations selected by the Bayesian Subset Simulation strategy at stage t = 1,2,3 and the final 
Design of Experiment (DoE). 



To evaluate the statistical properties of the estimator, we consider the absolute relative bias 

E(S) — a 



a 



and the coefficient of variation 



cov 



5(a) 



Q: 



(9) 



(10) 



where E(a) is the average mean and 5(a) is the standard deviation of the estimator a. 



Table 4 shows the results of the comparison of our proposed Bayesian Subset Sampling algorithm with the 
Subset Simulation algorithm in [3]. Crude Monte Carlo sampling is used as the reference probability of failure. 
For the fairness of the comparison, we set the same intermediate probability po = 0.1, and sample size m = 
1000 in both methods. Fifty independent runs are performed to evaluate the average of both methods. For 
Bayesian Subset Simulation method, as N is different for each run, we show the minimal and maximal of N 
from 50 runs. 



7 



Table 4: Compare with Monte Carlo approach and Subset Simulation 



Method 


m 


N 


a (10~ 5 ) 


5{a) (KT 5 ) 


K 


cov 


MCS 


10 8 


10 8 


3.8500 


0.062 





1.6% 


Sub-Sim 


1000 


4600 


3.9078 


2.470 


1.5% 


63.2% 


Bayesian Sub-Sim 


1000 


[94, 109] 


3.7020 


0.618 


4.4% 


16.7% 



4. CONCLUSION 

In this paper, we propose a new algorithm called Bayesian Subset Simulation for estimating small probabilities 
of failure in a context of very expensive simulations. This algorithm combines the main ideas of the Subset 
Simulation algorithm and the SUR strategies developed in our recent work [2]. 

Our preliminary results show that the number of evaluations is dramatically decreased compared to the original 
Subset Simulation algorithm, while keeping a small bias and coefficient of variation. 

Our future work will try to improve further the properties of our algorithm regarding the bias and the variance 
of the estimator. We shall also test and validate the approach on more challenging examples. 

ACKNOWLEDGMENTS 

The research of Ling Li, Julien Beet and Emmanuel Vazquez was partially funded by the French Fond Unique 
Interministeriel (FUI) in the context of the project CSDL. 

REFERENCES 

[1] E. Vazquez and J. Beet. A sequential Bayesian algorithm to estimate a probability of failure. In Proceed- 
ings of the 15th IFAC Symposium on System Identification, SYSID 2009 15th IFAC Symposium on System 
Identification, SYSID 2009, Saint-Malo France, 2009. 

[2] J. Beet, D. Ginsbourger, L. Li, V Picheny, and E. Vazquez. Sequential design of computer experiments 
for the estimation of a probability of failure. Statistics and Computing, pages 1-21, 2010. 

[3] S. K. Au and J. Beck. Estimation of small failure probabilities in high dimensions by subset simulation. 
Probab. Engrg. Mechan., 16(4):263-277 ', 2001. 

[4] F Cerou, P. Del Moral, T Furon, and A. Guyader. Sequential monte carlo for rare event estimation. 
Statistics and Computing, pages 1-14, 2011. 

[5] J. Sacks, W. J. Welch, T J. Mitchell, and H. P. Wynn. Design and analysis of computer experiments. 
Statistical Science, 4(4):409-435, 1989. 

[6] C. Currin, T Mitchell, M. Morris, and D. Ylvisaker. Bayesian prediction of deterministic functions, with 
applications to the design and analysis of computer experiments. J. Amer. Statist. Assoc., 86(416):953- 
963, 1991. 

[7] W. J. Welch, R. J. Buck, J. Sacks, H. P. Wynn, T J. Mitchell, and M. D. Moms. Screening, predicting and 
computer experiments. Technometrics, 34:15-25, 1992. 

[8] J. Oakley and A. O'Hagan. Bayesian inference for the uncertainty distribution of computer model outputs. 
Biometrika, 89(4), 2002. 



8 



[9] J.E. Oakley and A. O'Hagan. Probabilistic sensitivity analysis of complex models: a Bayesian approach. 
Journal of the Royal Statistical Society: Series B (Statistical Methodology), 66(3):75 1-769, 2004. 

[10] J. Oakley. Estimating percentiles of uncertain computer code outputs. /. Roy. Statist. Soc. Ser. C, 53(1):83- 
93, 2004. 

[11] M. J. Bayarri, J. O. Berger, R. Paulo, J. Sacks, J. A. Cafeo, J. Cavendish, C.-H. Lin, and J. Tu. A 
framework for validation of computer models. Technometrics, 49(2): 138-154, 2007. 

[12] J. Mockus, V. Tiesis, and A. Zilinskas. The application of Bayesian methods for seeking the extremum. 
In L. Dixon and Eds G. Szego, editors, Towards Global Optimization, volume 2, pages 1 17-129. Elsevier, 
1978. 

[13] J. Mockus. Bayesian Approach to Global Optimization. Theory and Applications. Kluwer Academic 
Publisher, Dordrecht, 1989. 

[14] D. R. Jones, M. Schonlau, and J. William. Efficient global optimization of expensive black-box functions. 
Journal of Global Optimization, 13(4):455^192, 1998. 

[15] J. Villemonteix. Optimisation defonctions couteuses. PhD thesis, Universite Paris-Sud XI, Faculte des 
Sciences d'Orsay, 2008. 

[16] J. Villemonteix, E. Vazquez, and E. Walter. An informational approach to the global optimization of 
expensive-to-evaluate functions. Journal of Global Optimization, 44(4): 509-534, 2009. 

[17] D. Ginsbourger. Metamodeles multiples pour V approximation et V optimisation de fonctions numeriques 
multivariables. PhD thesis, Ecole nationale superieure des Mines de Saint-Etienne, 2009. 

[18] V. Dubourg, B. Sudret, and J.-M. Bourinet. Reliability-based design optimization using kriging surrogates 
and subset simulation. Structural Multisciplinary Optimization, 44(5):673-690, 2011. 

[19] V. Dubourg. Adaptive surrogate models for reliability analysis and reliability-based design optimization. 
PhD thesis, Universite Blaise Pascal - Clermont II, 2011. 

[20] Pierre Del Moral, Arnaud Doucet, and Ajay Jasra. Sequential monte carlo samplers. Journal of the Royal 
Statistical Society: Series B (Statistical Methodology), 68(3):41 1^136, 2006. 

[21] R. Douc and E. Moulines. Limit theorems for weighted samples with applications to sequential monte 
carlo methods. The Annals of Statistics, 36(5):2344-2376, 2008. 

[22] Vincent Dubourg, Francois Deheeger, and Bruno Sudret. Metamodel-based importance sampling for 
structural reliability analysis. Preprint submitted to Probabilistic Engineering Mechanics, 201 1. 

[23] Vincent Dubourg, Francois Deheeger, and Bruno Sudret. Metamodel-based importance sampling for the 
simulation of rare events. In 11th International Conference on Applications of Statistics and Probability 
in Civil Engineering (1CASP 11), 2011. 

[24] Christian P. Robert and G. Casella. Monte Carlo statistical methods, 2nd edition. Springer Verlag, 2004. 

[25] T. J. Santner, B. J. Williams, and W. Notz. The Design and Analysis of Computer Experiments. Springer 
Verlag, 2003. 

[26] J. D. Hoi, T. B. Schon, and F. Gustafsson. On resampling algorithms for particle filters. In IEEE Workshop 
on Nonlinear Statistical Signal Processing Workshop, pages 79-82, 2006. 



9 



[27] Randal Douc and Olivier Cappe. Comparison of resampling schemes for particle filtering. In Proceedings 
of the 4th International Symposium on Image and Signal Processing and Analysis (ISPA), pages 64-69, 
2005. 

[28] M. Bolic, Petar M. Djuric, and S. Hong. New resampling algorithms for particle filters. In IEEE In- 
ternational Conference on Acoustics, Speech, and Signal Processing, 2003. Proceedings. (ICASSP '03), 
volume 2, pages 589-592, 2003. 

[29] N. Gayton, J.M. Bourinet, and M. Lemaire. Cq2rs: a new statistical approach to the response surface 
method for reliability analysis. Structural Safety, 25(1):99— 121, 2003. 

[30] M.R. Rajashekhar and B.R. Ellingwood. A new look at the response surface approach for reliability 
analysis. Structural Safety, 12(3):205-220, 1993. 

[31] Jason L. Loeppky, Jerome Sacks, and William J. Welch. Choosing the sample size of a computer experi- 
ment: A practical guide. Technometrics, 51(4):366-376, 2009. 

[32] M. L. Stein. Interpolation of Spatial Data: Some Theory for Kriging. Springer, New York, 1999. 



10 



